Chapter 7 



Spectral Interpolation, 
Differentiation, Quadrature 

7.1 Interpolation 

7.1.1 Bandlimited interpolation 

While equispaced points generally cause problems for polynomial interpola- 
tion, as we just saw, they are the natural choice for discretizing the Fourier 
transform. For data on xj = jh, j G Z, recall that the semidiscrete Fourier 
transform (SFT) and its inverse (ISFT) read 

/>) = hY^ e-'^^'fj, fj = ^ e^'^'fik) dk. 
^ 27r 7-./^ 

The idea of spectral interpolation, or bandlimited interpolation, is to evaluate 
the ISFT formula above at some point x not equal to one of the Xj. 

Definition 21. (Fourier/ spectml/handlimited interpolation on Let Xj = 
jh, j e Z. Consider / : R i— > its restriction fj = f{xj), and the SFT f{k) 
of the samples fj. Then the spectral interpolant is 

p{x) = ^ e^'^fik) dk. 

We can view the formula for p{x) as the inverse Fourier transform of the 
compactly supported function equal to f{k) for k e [— 7r//i, 7r//i], and zero 
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otherwise. When the Fourier transform of a function is compactly supported, 
we say that that function is bandlimited, hence the name of the interpolation 
scheme. 

Example 22. Let 

f - A - / 1 '^•^■ = 0' 

Then the SET is f{k) = h for k G [—n/h, n/h] as we saw previously. Extend 
it to k & M. by zero outside of [—■n/h, ■n/h]. Then 

s\n{'Kx/h) 
pix) — — = sincinx/h). 

TTX/h 

This function is also called the Dirichlet kernel. It vanishes at xj — jh for 
j 7^ 0; integer. 

Example 23. In full generality, consider now the sequence 

k&Z 

By linearity of the integral, 

Pi^) = X]/fcSznc(7r(x- - Xk)/h). 
fcez 

The interpolant is a superposition of sine functions, with the samples fj as 
weights. Here sine is the analogue of the Lagrange elementary polynomials 
of a previous section, and is called the interpolation kernel. For this reason, 
bandlimited interpolation sometimes goes by the name Fourier-sine interpo- 
lation. 

(Figure here for the interpolation of a discrete step.) 

In the example above we interpolate a discontinuous function, and the 
result is visually not very good. It suffers from the same Gibbs effect that 
we encountered earlier. The smoother the underlying f{x) which fj are the 
samples of, however, the more accurate bandlimited interpolation. 

In order to study the approximation error of l)audlimited interpolation, 
we need to return to the link between SFT and FT. The relationship be- 
tween p{x) and fj is sampling, whereas the relationship between the FT 
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f{k)x[—K/h,Tv/h]{k) and the SFT f{k) is periodization. We have already aUuded 
to this correspondence ear her, and it is time to formulate it more precisely. 

(Figure here; sampling and periodization) 

Theorem 14. (Poisson summation formula, FT version) Let -u : M i— >■ 
sufficiently smooth and decaying sufficiently fast at infinity. (We are delib- 
erately imprecise!) Let Vj — u{xj) for Xj — jh, j e and 

u{k)= ! e-'^''u{x)dx, {FT), A; e R, 
v{k) ^h^e-'^''^u{xj). (SFT), k e [-7r/h,7r/h]. 

Then 

v[k) = ^u{k + m—), k e [-tt/H, tt/H] (7.1) 

In some texts the Poisson summation formula is written as the special 
case A; = 0: 

h^u{jh) = ^«(-^m). 

jSZ meZ 

Exercise: use what we have already seen concerning translations and Fourier 
transforms to show that the above equation implies (hence is equivalent to) 
equation (7.1). 

Proof. Consider the right-hand side in (7.1), and call it 

0(A;) = ^ u[k + n^-^), k e [-7r//i, -k/H]. 

mez 

It suffices to show that ^{k) — v{k), or equivalently in terms of their ISFT, 
that (t)j = Vj for j e Z. The ISFT is written 
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The function u is smooth, hence intcgrable, and the sum over m converges 
fast. So we can interchange sum and integrah 

= y [ ^ U{k + TX?-^y^'^ dk. 

Now put k' — k + m^, and change variable: 

0, = — y / ' u{ky''''e-''^^''dk'. 



h h 

,-27r , 



The extra exponential factor e'~'^~^^^ is equal to 1 because j e Z. We are 
in presence of an integral over M chopped up into pieces corresponding to 
sub- intervals of length 27r/h. Piecing them back together, we get 



- / uiky^'^^dk', 



which is exactly the inverse FT of u evaluated at Xj = jh, i.e., (pj = u{xj) — 

Vj. 

The Poisson summation formula shows that sampling a function at rate 
h corresponds to periodizing its spectrum (Fourier transform.) with a period 
27r//i. So the error made in sampling a function (and subsequently doing 
bandlimited interpolation) is linked to the possible overlap of the Fourier 
transform upon periodization. 

• Scencirio 1. Assume supp(m) C [—n/h, n/h]. Then no error is made in 
sampling and interpolating u at rate h, because nothing happens upon 
27r//i-periodization and windowing into [—7r/h,7r/h]: 

p{k) = u{k) ^ 'p{x) = u{x). 

(Draw picture) 

• Scenario 2. Now assume that u{k) is not included in [—n/h, n/h]. 
In general the periodization of u{k) will result in some overlap inside 
[—n/hjir/h]. We call this aliasing. In that case, some information is 
lost and interpolation will not be exact. 
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Scenario 1 is known as the Shannon sampling theorem: a function ban- 
dhmited in [—7T/h,7T/h] in k space is perfectly interpolated by bandlimited 
interpolation, on a grid of spacing h or greater. In signal processing h is also 
called the sampling rate, because x has the interpretation of time. When h 
is the largest possible rate such that no ahasing occurs, it can be referred to 
as the Nyquist rate. 

More generally, the accuracy of interpolation is linked to the smoothness 
of u{x). If the tails u{k) are small and h is large, we can expect that the 
error due to periodization and overlap won't be too big. The following result 
is a consequence of the Poisson summation formula. 

Theorem 15. (Error of bandlimited interpolation) Let u have p > 1 deriva- 
tives in L^{M.). Let Vj — u{xj) at Xj — jh, j e Z. Denote by p{x) the 
bandlimited interpolant formed form Vj . Then, as h ^ 0, 

\u{k) - p{k)\ = 0{hP) \k\<^ 



and 

\\u-p\\2^0{hP-'/^). 



Proof. Denote by u{k) the FT of u{x), and by v{k) the SFT of vj, so that 
p{k) — v{k) on [— 7r//i, 7r//i]. By the Poisson summation formula (7.1), 

v{k) — u{k) = u{k + m—), k E [—n/h,Ti/h]. 

As we saw earlier, the smoothness condition on u imply that 

\u{k)\ <C\k\-P. 

Since 

27r TT 2-K TT 2tt, 
k + m— G — - + m— . - + m— , 
h h h h h 

we have \k + m'^\ > \mj^\, hence 

\u(k + m'^)\ < C'lm^r^ 
h h 

for some different constant C. Summing over m 7^ 0, 

■lyp < c"(- 



v{k)-u{k)\ <c'Y, H'^ily < c"C^)-p < C"'hP. 
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One can switch back to the x domain by means of the Plancherel formula 

ll^-^'lli^ = ^ll^(^) -'^(^)IIl2(r)- 

The right-hand side contains the integral of \u{k) — v{k)\'^ over R. Break this 
integral into two pieces: 

• Over [— 7r//i, 7r//i], we have seen that \u{k) — v{k)\'^ = 0{h'^P). The 
integral is over an interval of length 0(l/h), hence the norm squared 
is 0{h'^P~^). Taking a square root to get the norm, we get 

• For \k\ > n/h, we have p{k) = 0, so it suffices to bound J^^^^^f^ \'^ik)\'^ dk 
Since \u{k)\ < C \k\-P, this integral is bounded by 0{{7r/hfP-^). Tak- 
ing a square root, we again obtain a 

□ 

We have seen how to interpolate a function defined on R, but as a closing 
remark let us notice that a similar notion exists for functions defined on 
intervals, notably x e [— 7r,7r] or [0,27r]. In that case, wavenumbers are 
discrete, the FT is replaced by the FS, and the SFT is replaced by the DFT. 
Evaluating the DFT for x not on the grid would give an interpolant: 

N/2 
k=-N/2+l 

Contrast with the formula (6.8) for the IDFT. This is almost what we want, 
but not quite, because the highest wavenumber k — N/2 is treated asym- 
metrically. It gives rise to an unnatural complex term, even if fk is real and 
even. To fix this, it is customary to set f-N/2 = fN/2, to extend the sum 
from —N/2 to N/2, but to halve the terms corresponding to k — —N/2 and 
N/2. We denote this operation of halving the first and last term of a sum 
by a double prime after the sum symbol: 

E" 

Itis easy to check that this operation does not change the interpolating prop- 
erty. The definition of bandlimited interpolant becomes the following in the 
case of intervals. 
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Definition 22. (Spectral interpolation on [0, 2n]) Letxj = jh, j = 0, . . . , A^— 
1 with h = 1/N. Consider f : [0,2n] M., its restriction fj = f{xj), and 
the DFT fk of the samples fj. Then the spectral interpolant is 

p(x) = — > e'^-'fk. 

^ ' l-K^ k=-N/2 

Because p{x) is a superposition of "monomials" of the form e**^^ = (e^^)^ 
for k integer, we call it a trigonometric polynomial. 

The theory for interpolation by trigonometric polynomials (inside [0, 27r]) 
is very similar to that for general bandlimited interpolants. The only impor- 
tant modification is that [0, 2tt] is a periodized interval, so a function qualifies 
as smooth only if it connects smoothly across x = indentified with x = 27i 
by periodicity. The Poisson summation formula is still the central tool, and 
has a counterpart for Fourier series. 

Theorem 16. (Poisson summation formula, FS version) Letu : [0, 27r] i— > M, 
sufficiently smooth. Let vj — u{9j) for 9j = jh, j = 1, . . . ,N, h = 2n/N, 
and 

/•27r 

u^^ e-'^^u{e)de, {FS), keZ, 
Jo 

N 

Vk = hY,e-'''^u{9,). {DFT), k ^ . . . ,^ - I. 

Then 

, N N , ^ 

Vk = 2^Uk+mN, k = - — -1. (7.2) 

(Recall that N = ^ so this formula is completely analogous to (7.1).) 

For us, the important consequence is that if u has p derivatives in L^, 
over the periodized interval [0,27r], then the bandlimited interpolation error 
is a 0{hP) in the pointwisc sense in k space, and 0{hP~^^'^) in the L^ sense. 

For this result to hold it is important that u has p derivatives at the 
origin as well (identified by periodicity with 27r), i.e., the function is equally 
smooth as it straddles the point where the interval wraps around by pe- 
riodicity. Otherwise, if u{6) has discontinuities such as. u{27r~) ^ u{0'^), 
interpolation will suffer from the Gibbs effect. 
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7.1.2 Chebyshev interpolation 

Consider now smooth functions inside [—1,1] (for illustration), but not nec- 
essarily periodic. So the periodization f{x + 2j) may be discontinuous. 
Polynomial interpolation on equispaced points may fail because of the Runge 
phenomenon, and bandlimited interpolation will fail due to Gibbs's effect. 

A good strategy for interpolation is the same trick as the one we used for 
truncation in the last chapter: pass to the variable 9 such that 

X — cos 6. 

Then x e [—1, 1] corresponds to ^ e [0, tt]. We define g{9) = /(cos 6') with g 
27r-periodic and even. 

We can now consider the bandlimited interpolant of g{9) on an equispaced 
grid covering [0, 27r], like for instance 

TT? 

Oj^^, wherej = l,...,2iV. 

Using the definition we saw at the end of the last section (with the " double 
prime"), we get 

j=o 

By even symmetry in 9 and k (why?), we can write 

N 

^(^) ^ ^cos{k9)ck, 

k=0 

for some coefficients that are formed from the samples g{9j). 

Back to X, we get the sample points xj = cos{9j). They are called Cheby- 
shev points. They are not equispaced anymore, and bccausct they are the 
projection on the x-axis of equispaced points on the unit circle, they cluster 
near the edges of [— 1, 1] . There are N + 1 oi them, from j = to A^, because 
of the symmetry in 9. It turns out that the Xj are the extremal points of the 
Chebyshev polynomial Tn{x), i.e., the points in [—1,1] where T^ takes its 
maximum and minimum values. 
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In terms of the variable x, the interpolant can be expressed as 

N 

p{x) — gf(acos x) — ^^T„(x)c„, Tn{x) — cos(nacos x), 

n=0 

with the same c„ as above. Since T„(a;) are polynomials of degree < n < N, 
and p{x) interpolates f{x) at the + 1 (non-equispaced) points Xj, we are 
in presence of the Lagrange interpolation polynomial for / at Xj\ 

Now, the interesting conclusion is not so much the new formula for the in- 
terpolation polynomial in terms of T„, but that although this is a polynomial 
interpolant, the error analysis is inherited straight from Fourier analysis. If 
/ is smooth, then g is smooth and periodic, and we saw that the bandlimited 
interpolant converges very fast. The Chebyshev interpolant of / is equal to 
the bandlimited interplant of g so it converges at the exact same rate (in L°° 
in /c or n space) — for instance 0{N~p~^) when / has p derivatives (in BV). 

In particular, we completely bypassed the standard analysis of error of 
polynomial interpolation, and proved universal convergence for smooth /. 
The factor 



that was posing problems in the error estimate then does not pose a problem 
anymore, because of the very special choice of Chebyshev points cos{7i j/N) 
for the interpolation. Intuitively, clustering the grid points near the edges 
of the interval [—1,1] helps giving 7rN+i{x) more uniform values throughout 
[— 1, 1], hence reduces the gross errors near the edges. 

Let us now explain the differences in the behavior of the monic polynomial 
Y[f=oi^~^j) for equispaced vs. Chebyshev points, and argue that Chebyshev 
points are near-ideal for interpolation of smooth functions in intervals. The 
discussion below is mostly taken from Trefethen, p. 43. (See also the last 
problem on homework 2 for an example of different analysis.) 

Let p{z) — Y[j=oi^ ~ -^i)' where we have extended the definition of the 
monic polynomial to z e C. We compute 



TV 



j=0 



N 




j=0 
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Put 

N 

^n{z) = {N + l)-^^\og\z-Xj\. 

j=0 

The function 0Ar is like an electrostatic potential, due to charges at ^ = Xj, 
each with potential {N + l)~^log \ z — Xj\. Going back to p{z) from ^at is 
easy: 

Already from this formula, we can see that small variations in will lead 
to exponentially larger variations in p{z), particularly for large N. 

Let us now take a limit N ^ oo, and try to understand what happens 
without being too rigorous. What matters most about the Chebyshev points 
is their density: the Chebyshev points are the projection onto the real-axis 
of a sequence of equispaced points on a circle. If the density of points on the 
circle is a constant l/(27r), then the density of points on [—1, 1] generated by 
vertical projection is 

Pcheb{x) — — . (normalized to integrate to 1 on [—1,1]) 

TTvl — 

(This is a density in the sense that 

N / pcheb{x)dx 

J a 

approximately gives the number of points in [a, b].) Contrast with a uniform 
distribution of points, with density 

Pequiip) — 2' 

Then the potential corresponding to any given p{x) is simply 
(p{z) — j p{x)\og\z — x\dx. 

The integral can be solved explicitly for both densities introduced above: 
• For pequi, we get 

(t>equi{z) = -1 + ]^Re{{z + 1) \0g{z + l) - (z - 1) \og{z - 1)). 
It obeys 0egm(O) = -1, 0egm(±l) = "1 + log 2. 
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• For pcheb, we get 

I z — 's/ zP' — 1 
(l>Cheb{z) = log 

This function obeys (interesting exercise) (pchebix) = —log 2 for all 
X e [—1, 1] on the real axis. 

The level curves of both 4>equiv and (pcheb in the complex plane are shown 
on page 47 of Trefethen. 

Overlooking the fact that we have passed to a continuum limit for the 
potentials, we can give a precise estimate on the monic polynomial: 

(iv+i)^,,„i(^) _ / (2/e)^ near x = ±1; 



\Ve<im\Z)\ -e - j ^^^^^^ ^^^^ 

whereas 

\pcheb{z)\ ^ e(^+^)'^--(^) = 2-^, z e [-1, 1]. 

We see that Pequi can take on very different values well inside the interval 
vs. near the edges. On the other hand the values of pcheb are near-constant in 
[— 1, 1]. The density pchebi^) is the only one that will give rise to this behav- 
ior, so there is something special about it. It is the difference in asymptotic 
behavior of (2/e)~^ vs. that makes the whole difference for interpola- 
tion, as iV — > cxD. 

One may argue that neither Pequi nor pcheb blow up as A?" — > oo, but it 
is an interesting exercise to show that if we were interpolating in an interval 
[—a, a] instead of [—1,1], then the bounds would be multiplied by , by 
homogeneity. 

The Chebyshev points are not the only one that correspond to the density 
Pcheb{x) as A?" — > oo. For instance, there is also the Chebyshev roots 

which are the roots of T/v(,t), instead of being the extremal points. They 
give rise to very good interpolation properties as well. 

Finally, let us mention that the theory can be pushed further, and that 
the exponential rate of convergence of Chebyshev interpolation for analytic 
functions can be linked to the maximum value of (j){z) on the strip in which 
the extension f{z) of f{x) is analytic. We will not pursue this further. 
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7.2 Differentiation 

7.2.1 Bandlimited differentiation 

The idea that a numerical approximation of a derivative can be obtained from 
differentiating an interpolant can be pushed further. In this section we return 
to bandlimited interpolants. We've seen that they are extremely accurate 
when the function is smooth and periodic; so is the resulting differentiation 
scheme. It is called bandlimited differentiation, or spectral differentiation. 

First consider the case of a; G M and xj = jh, j £ Z. As we've seen, the 
bandlimited/spectral interpolant of u{xj) is 



where v{k) is the SFT of vj = u{xj). Differentiating p{x) reduces to a 
multiplication by ik in the Fourier domain. Evaluating p'{xj) is then just a 
matter of letting x = xj in the resulting formula. The sequence of steps for 
bandlimited differentiation [x e M) is the following: 

• Obtain the SFT v{k) of Vj — u{xj); 

• Multiply v{k) by ik; 

• Obtain the ISFT of w{k) — ikv{k), call it wj. 

The numbers Wj obtained above are an approximation of u'{xj). The 
following result makes this precise. 

Theorem 17. (Accuracy of bandlimited differentiation, see also Theorem 4 
in Trefethen's book) Let u have p derivatives in L-^(R). Let Vj — u{xj), and 
Wj —p'{xj) be the result of bandlimited differentiation. Then 




sup \wj — u'{xj) \ = 0{h^ ^). 



3 



and 



u — 



Proof. The proof hinges on the fact that 



v{k)-u{k)\ = 0{hP). 
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One power of h is lost when differentiating u (because ik is on the order of 
1/h over the fundamental cell [—Tr/h, 'n'/h]). Half a power of h is lost in going 
back to the physical domain {j instead of k) via the norm (why?), and a 
full power of h is lost when going back to j in the uniform sense (why?). □ 

The point of the above theorem is that the order of bandlimited differ- 
entiation is directly linked to the smoothness of the function, and can be 

arbitrarily large. This is called spectral accuracy. One can even push the 
analysis further and show that, when / is real-analytic, then the rate of 
convergence of Wj towards u'{xj) is in fact exponential/geometric. 

Of course in practice we never deal with a function u{x) defined on the 
real line. In order to formulate an algorithm, and not simply a sequence of 
abstract steps, we need to limit the interval over which u is considered. Spec- 
tral differentiation in the periodic interval [0, 27r] works like before, except 
DFT are substituted for SFT. For 9 e [0,27r], we've seen that the spectral 
interpolant is defined as 



(The double prime is important here.) Again, a derivative can be imple- 
mented by multiplying by ik in the Fourier domain. The sequence of steps 
is very similar to what it was before, except that we can now label them as 
"compute" , and not just "obtain" : 

• Compute the DFT Vk of Vj = u{xj); 

• Multiply Vk by ik; 

• Compute the IDFT of Wk — ikvk, call it Wj. 

The result of accuracy are the same as before, with the provision that u needs 
to be not only smooth, but also smooth when extended by periodicity. 

The FFT can be used to yield a fast 0{N log N) algorithm for spectral 
differentiation. 

Note that higher derivatives are obtained in the obvious manner, by mul- 
tiplying in Fourier by the adequate power of ik. 
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7.2.2 Chebyshev differentiation 

In view of what has been covered so far, the idea of Chebyshev differentiation 
is natural: it is simply differentiation of the Chebyshev interpolant at the 
Chebyshev nodes. It proves very useful for those smooth functions on an 
interval, which do not necessarily extend smoothly by periodicity 

Let us recall that the Chebyshev interpolant is q{x) = p(arccosa;), where 
p{9) is the bandlimited interpolant of u{cos9j) at the Chebyshev points 
Xj — cos9j. As such, a differentiation on q{x) is not the same thing as a 
differentiation on p{9). Instead, by the chain rule, 

q'{x) = p(arccosx). 
— x^ 

The algorithm is as follows. Start from the knowledge of u{xj) at Xj — 
cos9j, 9j = jh, h = tt/N, and j = 1, . . . , N. 

• Perform an even extension of u{xj) to obtain u{cos9j) for 9j = jh, 
h — tt/N, and j = —N + 1, . . . , A^. Now we have all the equispaced 
sample of the periodic function u{cos9) for 9j covering [0, 27r], not just 
[0,7r]. 

• Take the DFT of those samples, call it v^, 

• Multiply hy ik, 

• Take the IDFT of the result ikvk, 

• Multiply by — 1/ ^Jl — x^ to honor the chain rule. At the endpoints 

Xj = —1 or 1, take a special limit to obtain the proper values of 1) 
and See Trefethen's book for the correct values. 

This is a fast algorithm since we can use the FFT for the DFT and IDFT. 

Since we are only a change of variables away from Fourier analysis in a pe- 
riodic domain, the accuracy of Chebyshev differentiation is directly inherited 
from that of bandlimited differentiation. We also have spectral accuracy. 

Note that higher derivatives can be treated similarly, by applying the 
chain rule repeatedly. 

In practice, Chebyshev methods are particularly useful for boundary- 
value problems (we'll come back to this), when all samples of a function are 
to be determined at once, and when we have the freedom of choosing the 
sample points. 
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7.3 Integration 

7.3.1 Spectral integration 

To go beyond methods of finite order, the idea of spectral integration is to 
integrate a bandhmited interpolant. This strategy yields very high accuracy 
when the function is smooth and periodic. 

Consider a function u{9), 9 G [0,27r]. Its samples are Vj = u{9j) with 
9j = jh, h = 2ti/N, and j = 1, . . . , A^. Form the DFT: 

TV 

and the bandlimited interpolant, 

^ ' 27r^ k=-N/2 " 

Integrating p{9) gives the remarkably simple following result. 

/ p{9)d9^VQ^hY^u{9j). 

Jo ,=1 

We are back to the trapezoidal rule! (The endpoints are identified, = 9n-) 
While we have already encountered this quadrature rule earlier, it is now 
derived from Fourier analysis. So the plot thickens concerning its accuracy 
properties. 

Specifically, we have to compare Vq to Uq = Jq^ u{9) d9, where Uk are the 
Fourier series coefficients of u{9) . The relationship between vq and Uq is the 
Poisson summation formula (7.2): 

^ Ar 27r 

The most important term in this sum is for m = 0, and our task is again 
to control the other ones, for m 7^ 0. The resulting accuracy estimate is the 
following. 
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Theorem 18. Assume u has p derivatives in L-^[0,27r], where [0, 27r] is con- 
sidered a periodic interval. (So it matters that the function connects smoothly 
by periodicity.) Then 



/ u{9)d9-hY^u{9j) = 0{hP). 

7=1 



Proof. By the Poisson summation formula, in the notations of the preceding 
few paragraphs, 

Mo - = ^ UmN- 

We have aheady seen that the smoothness properties of u are such that 

\uk\ < C\k\-P. 

So we have 

\uo -vo\<c ^ 

which is a 0{hP). □ 

As a conclusion, the trapezoidal rule is spectrally accurate (error O^h'^'^^) 
for all p > when u G C°^), provided the function to be integrated in 
smooth and periodic. If the function is in an interval but is for instance 
discontinuous upon periodization, then we revert to the usual 0(/i^) rate. So 
the true reason for the trapezoidal rule generally being 0(/i^) and not 0{h°°) 
is only the presence of the boundaries! 

An important example of such periodic smooth function, is a regular 
C°°function multiplied by a C°° window that goes smoothly to zero at the 
endpoints of [a, b], like for instance a member of a partition of unity. (This 
plays an important role in some electromagnetism solvers, for computing the 
radar cross-section of scatterers.) 



7.3.2 Chebyshev integration 

Like for Chebyshev differentiation, we can warp 9 into x by the formula 
9 = arccosx, and treat smooth functions that are not necessarily periodic. 
Integrating a Chebyshev interpolant gives rise to Chebyshev integration, also 
called Clenshaw-Curtis quadrature. 
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7.3. INTEGRATION 



Assume that f{x) is given for a; G [—1,1], otherwise rescale the x variable. 
The Chebyshev interpolant is built from the knowledge of / at the Chebyshev 



nodes Xj = cos9j, and takes the form 



n>0 



We have seen that a„ are obtained by even extension of f{xj), followed by 
an FFT where the sine terms are dropped. As a result, 



We compute 



/ p{x) dx^y2 / Tn{x) dx. 

/I pTT 
T„{x)dx^ j cos{n0)sm0de. 



This integral can be evaluated easily by relating cos and sin to complex 
exponentials (or see Trefethen's book for an elegant shortcut via complex 
functions), and the result is 



J Tn{x) dx = 



if n is odd; 
if n is even. 



The algorithm for Chebyshev differentiation is simply: (1) find the coef- 
ficients an ioT < n < N, and (2) form the weighted sum 



E 



2 



"n ■ 



' -1 9 

1 — n'' 

n evcn,n<N 

The accuracy estimate of Chebyshev integration is the same as that of 
bandlimited integration, except we do not need periodicity. So the method 
is spectrally accurate. 

An important method related to Chebyshev quadrature is Gaussian quadra- 
ture, which is similar but somewhat different. It is also spectrally accurate for 
smooth functions. Instead of using extrema or zeros of Chebyshev polynomi- 
als, it uses zeros of Lcgcndrc polynomials. It is usually derived directly from 
properties of orthogonal polynomials and their zeros, rather than Fourier 
analysis. This topic is fascinating but would take us a httle too far. 
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